/**
 * @file   project.h
 * @author ubuntu <dadi@ubuntu>
 * @date   Sat Jun  5 22:33:51 2021
 * 
 * @brief  project2
 * 
 * 
 */

#ifndef __DADI__PROJECT2__
#define __DADI__PROJECT2__

#include<iostream>
#include<vector>
#include<tuple>
#include<map>
#include<string>
#include<cmath>

class Multigrid;

class Multigrid
{
protected:
    int n_grid;			/**< The number of grid.  */
    double h;			/**< The width of grid. */
    std::vector<double> RHS;	/**< The original right-hand side. */
    std::string reop;		/**< The restriction_operater pattern. */
    std::string inop;		/**< The interpolation_operater pattern. */
public:
    /** 
     * Set n_grid and h.
     * 
     * @param n 
     */
    void set_n_grid(int n);
    /** 
     * Set original RHS with the boundary condition.
     * 
     * @param boundary 
     * @param f 
     */
    void set_RHS(std::vector<double> boundary, double (*f)(double x));
    /** 
     * Get RHS.
     * 
     * 
     * @return The original RHS.
     */
    std::vector<double> get_RHS();
    /** 
     * The coefficient matrix operater.
     * 
     * @param u 
     * 
     * @return A*u.
     */
    std::vector<double> A_operater(std::vector<double> u);
    /** 
     * Compute the max_norm of a vector.
     * 
     * @param u 
     * 
     * @return 
     */
    double vector_max_norm(std::vector<double> u);
    /** 
     * Compute the discrete L2 norm of a vector.
     * 
     * @param u 
     * 
     * @return 
     */
    double discrete_L2_norm(std::vector<double> u);
    /** 
     * Weighted Jacobi method.
     *  
     * @param w weight coefficient
     * @param u 
     * @param f 
     * 
     * @return 
     */
    std::vector<double> weighted_Jacobi(double w, std::vector<double> u, std::vector<double> f);
    /** 
     * Map a vector from a fine grid(h) to a coarse grid(2h).
     * 
     * @param u 
     * 
     * @return 
     */
    std::vector<double> restriction_operater(std::vector<double> u);
    /** 
     * Map a vector from a coarse grid(2h) to a fine grid(h).
     * "full-weighting" and "injection"
     * @param u 
     * 
     * @return 
     */
    std::vector<double> interpolation_operater(std::vector<double> u);
    /** 
     * Pure virtual function, Multigridsolver.
     * "linear" and "quadratic"
     * @param init 
     * @param f 
     * @param w 
     * @param v1 
     * @param v2 
     * 
     * @return solution.
     */
    virtual std::vector<double> solver(std::vector<double> init, std::vector<double> f, double w, int v1, int v2) = 0; 
};

void Multigrid::set_n_grid(int n)
{
    n_grid = n;
    h = 1.0/n;
};

void Multigrid::set_RHS(std::vector<double> boundary, double (*f)(double x))
{
    RHS.resize(n_grid-1);
    for(int i = 0; i < n_grid-1; i++)
    {
	if(i == 0)
	    RHS[i] = f((i+1)*h) + boundary[0]/(h*h);
	else if(i == n_grid-2)
	    RHS[i] = f((i+1)*h) + boundary[1]/(h*h);
	else
	    RHS[i] = f((i+1)*h);
    }
};

std::vector<double> Multigrid::get_RHS()
{
    return RHS;
};

std::vector<double> Multigrid::A_operater(std::vector<double> u)
{
    int n = u.size() + 1;
    double _h;
    _h = 1.0/n;
    std::vector<double> u_new(n-1, 0.0);
    for(int i = 0; i < n-1; i++)
    {
	if(i == 0)
	    u_new[i] = (2.0*u[i] - 1.0*u[i+1])/(_h*_h);
	else if(i == n-2)
	    u_new[i] = (-1.0*u[i-1] + 2.0*u[i])/(_h*_h);
	else
	    u_new[i] = (-1.0*u[i-1] + 2.0*u[i] - 1.0*u[i+1])/(_h*_h);
    }
    u = u_new;
    return u;
};

double Multigrid::vector_max_norm(std::vector<double> u)
{
    double max_norm = 0;
    for(int i = 0; i < u.size(); i++)
	if(fabs(u[i]) > max_norm)
	    max_norm = fabs(u[i]);
    return max_norm;
};

double Multigrid::discrete_L2_norm(std::vector<double> u)
{
    double L2_norm = 0;
    int n = u.size() + 1;
    double h = 1.0/n;
    for(int i = 0; i < n-1; i++)
	L2_norm += h*u[i]*u[i];
    L2_norm = sqrt(L2_norm);
    return L2_norm;
};

std::vector<double> Multigrid::weighted_Jacobi(double w, std::vector<double> u, std::vector<double> f)
{
    int n = u.size() + 1;
    double _h = 1.0/n;
    std::vector<double> u_temp(n-1,0.0);
    u_temp = A_operater(u);
    for(int i = 0; i < n-1; i++)
	u_temp[i] = -w*_h*_h*u_temp[i]/2.0 + u[i] + w*_h*_h*f[i]/2.0;
    u = u_temp;
    return u;
};

std::vector<double> Multigrid::restriction_operater(std::vector<double> u)
{
    int n = u.size() + 1;
    std::vector<double> u_new(n/2 - 1, 0.0);
    if(reop == "full-weighting")
    {
	for(int i = 1; i < n/2; i++)
	    u_new[i-1] = (u[2*i-2] + 2.0*u[2*i-1] + u[2*i])/4.0;
    }
    else if(reop == "injection")
    {
	for(int i = 1; i < n/2; i++)
	    u_new[i-1] = u[2*i-1];
    }
    else;
    return u_new;
};

std::vector<double> Multigrid::interpolation_operater(std::vector<double> u)
{
    int n = 2*(u.size() + 1);
    std::vector<double> u_new(n-1, 0.0);
    if(inop == "linear")
    {
	for(int i = 1; i < n; i++)
	{
	    if(i%2 == 0)
		u_new[i-1] = u[i/2-1];
	    else if(i == 1)
		u_new[i-1] = u[(i-1)/2]/2.0;
	    else if(i == n-1)
		u_new[i-1] = u[(i-1)/2-1]/2.0;
	    else
		u_new[i-1] = (u[(i-1)/2-1] + u[(i-1)/2])/2.0;
	}
    }
    else if(inop == "quadratic")//Langrange interpolation
    {
	for(int i = 1; i <= n/4; i++)
	{
	    if(i == 1)
	    {
		u_new[4*i-4] = 3.0*u[2*i-2]/4.0 - u[2*i-1]/8.0;
		u_new[4*i-3] = u[2*i-2];
		u_new[4*i-2] = 3.0*u[2*i-2]/4.0 + 3.0*u[2*i-1]/8.0;
		u_new[4*i-1] = u[2*i-1];
	    }
	    else if(i == n/4)
	    {
		u_new[4*i-5] = u[2*i-3];
		u_new[4*i-4] = 3.0*u[2*i-3]/8.0 + 3.0*u[2*i-2]/4.0;
		u_new[4*i-3] = u[2*i-2];
		u_new[4*i-2] = - u[2*i-3]/8.0 + 3.0*u[2*i-2]/4.0;
	    }
	    else
	    {
		u_new[4*i-5] = u[2*i-3];
		u_new[4*i-4] = 3.0*u[2*i-3]/8.0 + 3.0*u[2*i-2]/4.0 - u[2*i-1]/8.0;
		u_new[4*i-3] = u[2*i-2];
		u_new[4*i-2] = - u[2*i-3]/8.0 + 3.0*u[2*i-2]/4.0 + 3.0*u[2*i-1]/8.0;
		u_new[4*i-1] = u[2*i-1];
	    }
	}
    }
    else;
    return u_new;
};

class V_cycle : public Multigrid
{
private:
public:
    /** 
     * Constructor.
     * 
     * @param _reop restriction_operater 
     * @param _inop interpolation_operater
     * 
     * @return 
     */
    V_cycle(std::string _reop, std::string _inop);
    /** 
     * V-cycle solver.
     * 
     * @param init 
     * @param f 
     * @param w 
     * @param v1 
     * @param v2 
     * 
     * @return 
     */
    std::vector<double> solver(std::vector<double> init, std::vector<double> f, double w, int v1, int v2) override;
};

V_cycle::V_cycle(std::string _reop, std::string _inop)
{
    n_grid = 0;
    h = 0;
    RHS = {0.0};
    reop = _reop;
    inop = _inop;
};

std::vector<double> V_cycle::solver(std::vector<double> init, std::vector<double> f, double w, int v1, int v2)
{
    int n = init.size() + 1;
    std::vector<double> u;
    std::vector<double> fh;
    u = init;
    fh = f;
    for(int i = 0; i < v1; i++)
	u = weighted_Jacobi(w, u, fh);
    if(n == 8)
    {
	for(int i = 0; i < v2; i++)
	    u = weighted_Jacobi(w, u, fh);
    }
    else
    {
	std::vector<double> f2h(n/2 - 1, 0.0);
	std::vector<double> v(n/2 - 1, 0.0);
	std::vector<double> temp(n-1, 0.0);
	for(int i = 0; i < n-1; i++)
	    temp[i] = fh[i] - A_operater(u)[i];
	f2h = restriction_operater(temp);
	v = solver(v, f2h, w, v1, v2);
	for(int i = 0; i < n-1; i++)
	    u[i] = u[i] + interpolation_operater(v)[i];
	for(int i = 0; i < v2; i++)
	    u = weighted_Jacobi(w, u, fh);
    }
    return u;
};

class Full_multigrid_V_cycle : public Multigrid
{
private:
public:
    /** 
     * Constructor.
     * 
     * @param _reop 
     * @param _inop 
     * 
     * @return 
     */
    Full_multigrid_V_cycle(std::string _reop, std::string _inop);
    /** 
     * Full multigrid V-cycle solver. Not perfect, init is not needed, 
     * just for the consistency of pure virtual function in the basic class.
     * @param init 
     * @param f 
     * @param w 
     * @param v1 
     * @param v2 
     * 
     * @return 
     */
    std::vector<double> solver(std::vector<double> init, std::vector<double> f, double w, int v1, int v2) override;
};

Full_multigrid_V_cycle::Full_multigrid_V_cycle(std::string _reop, std::string _inop)
{
    n_grid = 0;
    h = 0;
    RHS = {0.0};
    reop = _reop;
    inop = _inop;
};

std::vector<double> Full_multigrid_V_cycle::solver(std::vector<double> init, std::vector<double> f, double w, int v1, int v2)
{
    int n = init.size() + 1;
    std::vector<double> u(n-1, 0.0);
    std::vector<double> fh(n-1, 0.0);
    fh = f;
    if(n == 8)
    {
	V_cycle V(reop, inop);
	u = V.solver(u, fh, w, v1, v2);
    }
    else
    {
	std::vector<double> f2h(n/2 - 1, 0.0);
	std::vector<double> v(n/2 - 1, 0.0);
	std::vector<double> b(n/2 - 1, 0.0);
	f2h = restriction_operater(fh);
	v = solver(b, f2h, w, v1, v2);
	u = interpolation_operater(v);
	V_cycle V(reop, inop);
	u = V.solver(u, fh, w, v1, v2);
    }
    return u;
};

class MultigridFactory
{
public:
    typedef Multigrid* (*MultigridCreator)();

public:
    /** 
     * Register Multigrid pattern.
     * 
     * @param string 
     * @param string 
     * @param schema 
     * @param createFn 
     * 
     * @return 
     */
    bool RegisterMultigrid(const std::tuple<std::string, std::string, std::string> schema, MultigridCreator createFn);
    /** 
     * Unregister Multigrid pattern.
     * 
     * @param string 
     * @param string 
     * @param schema 
     * 
     * @return 
     */
    bool UnregisterMultigrid(const std::tuple<std::string, std::string, std::string> schema);
    /** 
     * Create a Multigrid pointer.
     * 
     * @param string 
     * @param string 
     * @param schema 
     * 
     * @return 
     */
    Multigrid* CreateMultigrid(const std::tuple<std::string, std::string, std::string> schema) const;

private:
    typedef std::map<std::tuple<std::string, std::string, std::string>, MultigridCreator> CallbackMap;
    CallbackMap callbacks_;

};

bool MultigridFactory::RegisterMultigrid(const std::tuple<std::string, std::string, std::string> schema, MultigridCreator createFn)
{
    return callbacks_.insert(CallbackMap::value_type(schema, createFn)).second;
};

bool MultigridFactory::UnregisterMultigrid(const std::tuple<std::string, std::string, std::string> schema)
{
    return callbacks_.erase(schema) == 1;
};

Multigrid* MultigridFactory::CreateMultigrid(const std::tuple<std::string, std::string, std::string> schema) const
{
// CallbackMap::const_iterator i = callbacks_.find(schema);
    auto i = callbacks_.find(schema);
    if (i == callbacks_.cend())
    { // not found
	throw std::runtime_error("Unknown schema");
    }
// Invoke the creation function
    return (i->second)();
};

Multigrid* CreateVC_fl()
{
    return new V_cycle("full-weighting", "linear");
};

Multigrid* CreateVC_fq()
{
    return new V_cycle("full-weighting", "quadratic");
};

Multigrid* CreateVC_il()
{
    return new V_cycle("injection", "linear");
};

Multigrid* CreateVC_iq()
{
    return new V_cycle("injection", "quadratic");
};

Multigrid* CreateFMG_fl()
{
    return new Full_multigrid_V_cycle("full-weighting", "linear");
};

Multigrid* CreateFMG_fq()
{
    return new Full_multigrid_V_cycle("full-weighting", "quadratic");
};

Multigrid* CreateFMG_il()
{
    return new Full_multigrid_V_cycle("injection", "linear");
};

Multigrid* CreateFMG_iq()
{
    return new Full_multigrid_V_cycle("injection", "quadratic");
};

#else
//DO NOTHING.
#endif
